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Motivated by recent experiments on atomic Dirac fermions in a tunable honeycomb optical lattice, 
we study the attractive Hubbard model of superfluidity in the anisotropic honeycomb lattice. At 
weak-coupling, we find that the maximum mean field pairing transition temperature, as a function 
of density and interaction strength, occurs for the case with isotropic hopping amplitudes. In this 
isotropic case, we go beyond mean field theory and study collective fluctuations, treating both 
pairing and density fluctuations for interaction strengths ranging from weak to strong coupling. 
We find evidence for a sharp sound mode, together with a well-defined Leggett mode over a wide 
region of the phase diagram. We also calculate the superfluid order parameter and collective modes 
in the presence of nonzero superfluid flow. The flow-induced softening of these collective modes 
leads to dynamical instabilities involving stripe-like density modulations as well as a Leggett-mode 
instability associated with the natural sublattice symmetry breaking charge-ordered state on the 
honeycomb lattice. The latter provides a non-trivial test for the experimental realization of the 
one-band Hubbard model. We delineate regimes of the phase diagram where the critical current is 
limited by depairing or by such collective instabilities, and discuss experimental implications of our 
results. 

PACS numbers: 03.75.Ss,71.10.Fd,81.05.ue,74.70.Wz 

I. INTRODUCTION 

Ultracold atoms in optical lattices are of great interest for studying, strongly correlated states of quantum matter, 
and such emergent phenomena as superconductivity and magnetism Q. Stimulated by the discovery of graphene, 
and ongoing intense research in that field [3|, optical lattices with the honeycomb lattice structure are beginning to 
be explored by various groups [4-6]. In the context of graphene, various exotic phenomena related to massless Dirac 
fermions as well as interesting topological phases and quantum phase transitions have been widely explored 
While superconductivity still remains elusive in experiments on graphene, it has been the focus of recent theoretical 
work suggesting chiral superconductivity induced by repulsive interactions [Iol-[l3|. Earlier, motivated by a possible 
cold atom realization, where the interatomic interaction may be modelled by an s-wave attractive contact potential, 
the attractive fermion Hubbard model on the two-dimensional (2D) honeycomb lattice was studied and found to 
exhibit a quantum phase transition at half-filling between a semimetal with massless Dirac fermion excitations and 
a simple s-wave superfluid phase [3, [H|. This quantum phase transition was shown to be the end point (with 
changing density) of a suitably defined BCS-BEC crossover line away from half filling [14[. This is reminiscent of 
the manner in which the usual BCS-BEC crossover phenomenon in the continuum in 3D at finite density is linked 
to the "zero density quantum critical point" associated with two-body bound state formation [16j. The semimetal 
to superconductor transition found on the honeycomb lattice is the simplest version of more general band insulator 
to superfluid transitions discussed in experiments [17j and theo ry |18l -[21|. The critical theory of this semimetal to 
superconductor transition is also a topic of great recent interest [22j . 

In this paper, motivated by the recent experimental developments reported by the ETH group [6], we present a 
careful study of Cooper pairing, superfluid collective modes, and superflow instabilities which limit the critical current 
in the attractive Hubbard model on the honeycomb lattice. Our main results, which go beyond previous work on 
this topic, are as follows, (i) We study the mean field pairing transition temperature T c ° in a general anisotropic 
honeycomb lattice as has been realized in a "brick- wall" type of geometry in the ETH experiments ^ , with hoppings 
as shown in Fig. [TJ For weak coupling, the maximum T c °, as a function of fermion density and interaction strength, 
is found to occur in the isotropic limit. We then argue more generally that the vicinity of this isotropic point, with 
all hoppings being equal, is the most promising limit to experimentally explore the superfluidity of atomic fermions. 
(ii) In this limit of isotropic hopping, we compute the collective modes in the superfluid, treating both density and 
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phase fluctuations on equal footing in contrast to previous work which only studied superconducting fluctuations [14[ . 
We find evidence for the expected sharp Anderson-Bogoliubov (AB) sound mode; in addition, over a wide swath of 
the phase diagram, we uncover a well-defined Leggett mode, which involves intra- unit-cell coupled oscillations of the 
density and superfluid order parameter phase, (iii) Imposing a nonzero superflow, we find that the superfluid phase 
becomes unstable via either a depairing instability or collective charge modulational instabilities driven by a softening 
of the collective mode spectrum. The superfluid is most stable against depairing in the vicinity of special fillings 
associated with van Hove (vH) singularities in the density of states; in this regime, at mean field level, we find the 
emergence of a supercurrent induced gapless superfluid state. The collective charge modulational instabilities that 
we uncover include tendencies towards forming incommensurate charge orders, stripe-like modulations, or sublattice 
symmetry breaking charge density wave (CDW) order. Aside from the depairing instability, which occurs primarily at 
small fermion density or in the weak attraction limit, we expect such collective instabilities to limit the critical current 
in this system. In fact, the current-induced gapless superfluid state we uncover in mean field theory is preempted by 
such collective dynamical or Landau instabilities for the intermediate coupling regime. Furthermore, the sublattice 
CDW order or the instability associated with this CDW order arises from a subtle pseudospin SU(2) symmetry of 
the attractive Hubbard model at half filling; probing for such CDW order can thus serve as a strong diagnostic tool 
to ascertain that the attractive one band Hubbard model is an appropriate description of the low energy physics in 
a deep optical lattice. Our predictions can be tested by using Bragg scattering experiments to study the collective 
mode spectrum. 

The results described above are obtained using a generalized random phase approximation (GRPA) complemented 
by strong coupling "pseudospin" -wave theory where appropriate. Our work here builds on similar previous studies 
addressing the collective modes and critical current-limiting instabilities on square and cubic optical lattices [23l-[25j. 
In those cases, it has been shown that the collective mode spectrum of s-wave superfluids exhibits, in addition to the 
long wavelength AB sound mode, a sharp 'roton minimum' near the corner of the Brillouin zone (BZ) due to strong 
CDW fluctuations, and 'roton mode softening' leads to a collective dynamical instability of superflow. 

This paper is organized as follows. In Sec. [Hj we first examine the honeycomb optical lattice setup from Ref. 
which can tune the anisotropy of Dirac cones, and discuss the mean field T c as a function of anisotropy A, interaction 
U, and filling n. In Sec. HIH we calculate the superfluid order parameter and quasiparticle dispersion in the presence 
of superflow. Section [IV] uses the GRPA to compute the spectrum of collective excitations as a function of flow. In 
the strong coupling limit, the GRPA result is compared with a spin wave expansion of the appropriate pseudospin 
Hamiltonian. GRPA works well even at strong coupling, as it shows excellent agreement with spin wave theory. 
Finally, Sec. IIVBI discusses the different kinds of superflow instabilities as a function of density and interaction 
strength. We find a dynamical instability at the T point which arises from competing CDW order. In addition, we 
find an unexpected stripe-like instability which persists even in the strong coupling limit. We conclude by discussing 
implications for experiments. 

II. HONEYCOMB LATTICE WITH TUNABLE ANISOTROPY 

The experiments by Tarruell et al. [6|, realize a tunable optical lattice which could be viewed as a "brick- wall" 
lattice or an anisotropic honeycomb lattice. In a deep optical lattice, we can restrict our attention to the lowest band 
and work with a single band tight-binding model. On symmetry grounds, two of the three nearest neighbor hoppings 
on this lattice are equivalent in the experiments, so we denote the three hopping amplitudes by (£, £, \t) as shown in 
Fig. [D (a). In the experiments, one pair of sites can be tuned to be closer to each other (A > 1) or further from each 
other (0 < A < 1); we therefore carry out a mean field study for general anisotropy A > 0. 

A. Attractive Hubbard model 

Since future experiments are likely to be able to study atomic fermions with attractive interactions in this lattice, 
we consider the attractive Hubbard model on this anisotropic honeycomb lattice as a reasonable starting point to 
study the superfluid phase. The Hamiltonian is thus given by 

H = ~ l V ( C i* C 3<r + C ]* C ™) ~ » U ™ ~ U n ^ n ^ ' w 

(i,j),a i,cr i 

where Cj a is the creation operator of a fermion with spin <r(=t, I) at site j, correspond to nearest neighbor sites 
coupled by a hopping amplitude t^-, the chemical potential fi tunes the fermion density, and U is the local Hubbard 
attraction. The hopping Uj = (t, t, Xt) is chosen as shown in Fig. Q] (a) for the three neighboring sites. 
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(a) 




FIG. 1: (Color online) (a) Honeycomb lattice in a "brick- wall" type of geometry. Sites of the A(B) sublattice are shown in 
green(red). Asymmetric hopping amplitudes are captured by a parameter A. A = 1 is the isotropic limit. The shaded area 
shows the unit cell, (b) Honeycomb lattice with primitive lattice vectors, (c) The Brillouin zone with the high symmetry points 
indicated. 



B. Single-particle dispersion 

The single-particle Hamiltonian in momentum space takes the form 

H = Y1 (^ C l,a,* C kA* + lt C l,b,a C k,a,a) (2) 

where the subscripts a, b refer to the two sublattices, and 

Ik = ~t(\ + e ika * + e ik < ai+a ^), (3) 

where a x = (y/3d,0) and a 2 = (- V^d/2, 3d/2) are the two basis vectors shown in Fig. [T] (b). This results in two 
dispersing bands with energies =b|7fc|. 

For U = 0, we plot the noninteracting band dispersion \jk\ for three different values of A as shown in Fig. [2j 
In the isotropic limit, A = 1, there are two zero energy Dirac cones which occur at the K points (BZ corners) as is 
familiar from graphene. Upon decreasing A from unity, the two Dirac points move towards the T point. When A <C 1, 
the system becomes quasi one dimensional and the Dirac points acquire highly anisotropic velocities, with a smaller 
velocity along the y direction. For A > 1, increasing A moves the Dirac cones along the BZ edge towards the M point 
(edge center), with a concomitant tunable velocity anisotropy. Eventually, at A = 2, the two Dirac points meet at the 
M point and 'annihilate' each other, resulting in a gapped spectrum [261428]. Such moving and merging/gapping of 
Dirac points have been observed in the ETH experiments 

Figure [3] shows the density of states (DOS) at the Fermi level as a function of A and fermion filling n (here n = 1 
corresponds to 'half filling' with one fermion on average per lattice site). In the isotropic A = 1 limit, the DOS 
exhibits vH singularities at n = 3/4,5/4 arising from the hexagonal shape of the Fermi surface at these densities. 




FIG. 2: (Color online) Top view of dispersion of the non- interacting fermions with A = 0.4 (left), A = 1 (center), and A = 1.6 
(right). The white dashed line is the boundary of the first Brillouin zone. In the symmetric A = 1 limit, the Dirac cones lie at 
the K points. 
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FIG. 3: (Color online) Left: g(cF), density of states at the Fermi level as a function of filling and A. We only show n < 1 as 
the system is particle- hole symmetric. Center: Mean field T<? /t with U/t = 1 as a function of filling (in units of fermions per 
site) and A. Right: Mean field T c °/t with U/t = 3 as a function of filling (in units of fermions per site) and A. 



Upon decreasing A from unity, we find that these vH peaks split into two branches. As we approach the A = limit, 
one of the branches merges with the 'band edge' at n = or n = 2. The other branch approaches n = 1, however it 
gets weaker and eventually vanishes at A = 0. Upon increasing A from unity, we again find that the vH singularities at 
n = 3/4, 5/4 split into two branches. For A ^> 1, the DOS peaks at densities n = 0, 2, at energies ±A£ corresponding 
to the bonding and antibonding state energies of isolated strong bonds. 



C. Cooper pairing: Mean field theory 

We next consider the effect of the attractive Hubbard interaction. At mean field level, the ground state of this 
model is generically a superfluid state formed by condensing Cooper pairs. To study this superfluid state, the fermions 
have to be cooled (at least) below the mean field transition temperature, T®, corresponding to the onset of Cooper 
pairing. We therefore begin by examining as a function of anisotropy A, and filling. 

At mean field level, the gap and number equations are given by [14[ 

1 i ^tanh^™ 



1 tanh^ 
U N ^ ^ 2E T (k) ' V ; 

k t=± 

and 
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respectively, where £^ = =b|7fc| — \i and E ± (k) = y ) 2 + Aq (r denotes the energy bands). N = 2M is the number 
of lattice sites and M is the number of unit cells. The derivation of Eqs. dlj and (|5]) is summarized in Appendix 
A. To determine T® (and the appropriate /i), we set Ao and /? -4- 1/T® in these equations and solve them 
self-consistently. 

Within weak-coupling BCS theory, the pairing gap and T c ° are determined by ^(e^), which is the DOS at the Fermi 
level in the non-interacting problem. We therefore compare the behavior of the mean field T® with the trends in 
g(ep)- At weak coupling (U/t = 1), we find that the maximum T c ° tracks g(ep) for all values of lambda and filling as 
seen from Fig. [3j As seen in the figure, the strong vH singularity in the symmetric A = 1 limit, leads to particularly 
robust in this isotropic case for n = 3/4,5/4. (We will revisit this robustness of the superfluidity in Sec lIII Bl 
in the context of imposed superflow.) For intermediate and large U/t, the pairing is no longer directly restricted to 
electrons near the Fermi surface, and the T c ° values are no longer governed simply by g{ep) as shown in Fig. fright). 
We find that for all interaction strengths, the mean field is highest when A < 1 as seen from Fig. [3j However, as 
the system becomes increasingly quasi- ID with decreasing A, we expect enhanced quantum and thermal fluctuations 
beyond mean field theory to suppress superfluidity, especially in the strong coupling regime where the superfluid to 
normal transition temperature is essentially controlled by fluctuation effects, and not by pair breaking effects which 
set T®. On these grounds, we expect the highest superfluid transition temperature to occur in the vicinity of the 
isotropic limit A = 1. We suggest A = 1 and n = 3/4 (or 5/4) as the optimal parameter regime for studying superfluid 
order in the experiments. In the rest of this paper, we therefore set A = 1 and study the collective excitations and 
superflow instabilities of the resulting superfluid state. 



III. NONZERO SUPERFLOW: MEAN FIELD THEORY ON THE ISOTROPIC HONEYCOMB LATTICE 



We next turn to the mean field theory of the isotropic honeycomb lattice Hubbard model in the presence of a 
nonzero superfluid current obtained by considering Cooper pairs with nonzero center of mass momentum Q. This 
will allow us to comment on the stability of the superfluid towards depairing (pair breaking) effects. This will also 
set the stage for the analysis of collective modes and collective mode driven instabilities discussed in later sections. 



A. Mean field theory in the presence of superfluid flow 



The order parameter in the presence of superfluid flow with flow momentum Q is given by 



&Q,»e iQri = Uid^ct^) 
U_ 
M 



J7 ^( C -k+Q/2,v,lCk+Q/2,v,t) elQ ' ri > ( 6 ) 



where I denotes the unit cell and v(= a, b) labels the sublattice. Within the mean field approximation, the Hamiltonian 
takes the form 



u 



where 



C k-,a,l 
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-/i 


- A Q,a 


7fc+ 


^ 




— A* 
















-fi 




V 





-7k- 


— A* 





(8) 



Here, k± denotes ±fc + Q/2. The matrix Hq can be diagonalized by the Bogoliubov transformation 



y Q (k)=f^ Q (k) = f 



(9) 
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where a\, T a is the creation operator of a Bogoliubov quasiparticle with momentum k + Q/2, energy band r(= ±), 
and spin <j(=t,i). They satisfy the anticommutation relation {a^^^a^, r , a ,} = ^fe'^T'^o-'- The unitary matrix 
T is determined so as to diagonalize the mean field Hamiltonian 



ft Q = fih Q f 



(E+{k) 
-E+(-k) 
££(fc) 





(10) 



(11) 



where Eq(k) and Eq(k) are the upper and lower quasiparticle energy bands, respectively. Thus, we obtain the 
diagonalized Hamiltonian 



k T = ± 



and the ground state energy 



E 



M 



#1' 



QM 



k T = ± 



(12) 



(13) 



To evaluate Aq a / fe and /i, we solve the corresponding gap and number equations that are explicitly given in Appendix 
lAl The chemical potential is tuned to obtain the required density and the order parameters are chosen to minimize the 
ground state energy. This scheme is known to interpolate between the weak-coupling BCS regime and strong-coupling 
BEC regime at low temperatures [29-31]. The strong-coupling effect is included in the self-consistently determined 
chemical potential, which reduces to the Fermi energy in the weak-coupling BCS regime while it deviates from the 
Fermi energy in the strong coupling BEC regime. 

For small order parameter and very small Q, we find a current in mean field theory which is linear in Q at all 
densities, including at and near half filling where a linearized Dirac theory [15, 32] would be applicable. However, our 
reported results mostly explore the regimes of large superflow momenta and a wide range of fillings, where a linearized 
Dirac dispersion is not applicable and the mean field renormalization of the gap must also be taken into account. 
Furthermore, as we show in the next section, once we consider physics beyond mean field theory, any nonzero current 
is unstable at half- filling for the attractive Hubbard model due to a dynamical instability associated with charge order. 



B. Depairing instability and gapless superfluidity 

Figured] shows the self-consistently calculated order parameter as a function of superflow momentum Q for different 
fillings. Hereafter in this section and in Sec. IIV| we restrict ourselves at T = for simplicity. For any flow direction, 
we find that the amplitudes of the two order parameters are the same, i.e., |Aq ? a| = |Aq,b| = Aq. We plot this 
amplitude Aq as a function of Q for different fillings in FigJU Since the Hamiltonian is particle-hole symmetric, 
we restrict ourselves to densities below half- filling (n < 1). Sweeping Q in the BZ, we find that the order parameter 
generally vanishes discontinuously at some critical flow momentum (except in the vicinity of n = 3/4). This depairing 
instability occurs when the kinetic energy cost of superflow outweighs the energy gain from pairing. The system then 
reverts to an unpaired normal state via a first order phase transition. 

However, when the filling is close to n = 3/4 fermions per site, Aq does not vanish for any Q. This enhanced 
stability is due to a vH singularity in the non-interacting problem at n = 3/4 and n = 5/4 (see Fig|3]). At these 
densities, there is a very large density of states at the Fermi level which enhances superfluid order and makes it robust 
against superflow. 

In Fig. 31 each curve for Aq exhibits an interesting cusp structure, at which the derivative OAq/OQ changes 
discontinuously. The order parameter starts decreasing rapidly beyond the cusp until the transition to the normal 
state sets in. This cusp structure signals the onset of gapless superfluidity. The evolution of the single-particle 
spectrum in the lower band Eq(k) is shown in Figs. [5] (a)~(c). With increasing flow, the energy gap to creating 
quasiparticles decreases for wavevectors k in the direction opposite to Q due to a 'Doppler shift' in the dispersion. 
For large Q, the gap closes making the single-particle spectrum gapless - this precisely coincides with the cusp feature. 
However, the closing of the gap does not immediately lead to depairing, as the self-consistent order parameter remains 
non-zero. 
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FIG. 4: (Color online) Superfluid order parameter Aq(= | Aq, | = | Aq^ |) as a function of the superflow Q for several different 
fermion fillings at T = for U/t = 4. Q is swept along the contour displayed in the inset. 
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FIG. 5: (Color online) ((a)~(c)) Single-particle spectrum E^fe) (in units of t) and ((d)~(f)) order parameter weight Wq A {k) 
(in units of t) in the lower excitation energy band for different flow momenta Q. We set U/t = 4, T = 0, and n = 0.5. In panel 
(c), Eq(k) takes negative values in the white regions indicating the gapless superfluid state. 



To clarify the origin of this cusp structure, we plot the single-particle spectrum of the lower energy band Eq(k), 
as well as the order parameter weight for the same band defined as 



= T* 3 T 13 f(EQ(k))+T* 4 T 14 (l - f(EQ(-k))). 



(14) 



Here, f(E) = \/{eP E + 1) is the Fermi distribution function and Tij denotes the ij th element of the Bogoliubov 
transformation matrix in Eq. (j9j). This quantity Wq a (k) is the contribution from wavevector k in the gap equation, 
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Eq. (|A1|) . As the quasiparticles are typically gapped, f{Eq(k)) = at T = and the only contribution to the weight 
comes from the term proportional to (1 — f(Eq(—k))). If the quasiparticle energy becomes negative for a certain 
momentum p, i.e., Eq(p) < 0, we have f(Eq(p)) = 1. Now, the order parameter weight at k = —p, as defined in 

Eq. (fl4|) . vanishes as (1 — f(Eq(—k))) = 0. This is clearly seen in Fig. |5] (c) and (f). The single-particle spectrum 
becomes gap less at some wavevectors p, and the order parameter weight vanishes at —p. Thus, the closing of the 
single-particle gap suppresses the order parameter weight at certain pockets. Consequently, the order parameter 
subsequently decreases sharply, but remains non-zero. This is the origin of the cusp in A versus Q at the onset of 
gapless superfluidity. 

IV. COLLECTIVE MODE ANALYSIS 

The mean field theory in Sec. IIII Al describes a superfluid state where Cooper pairs condense in a single momentum 
state with uniform density. To investigate the stability of superfluid order with imposed flow, we study quantum 
fluctuations beyond mean field theory. These fluctuations give rise to dispersing collective modes in the superfluid. 
For concreteness, we restrict ourselves to flow in the V — K direction in the rest of this paper. This direction is best 
suited for the optical lattice setup of Ref. 3 which uses three lasers at angles of 120° with respect to each other. Flow 
can easily be induced in the T — K direction by detuning one of the lattice lasers. For concreteness, we take flow to 
be parallel to the a\ lattice vector (see FigJTt)). Together with our definition of the unit cell, this forces Aq^ = Aq^ 
by symmetry. 

In general, the collective mode spectrum can be evaluated from the poles of the relevant response functions. In the 
superfluid state, we have fluctuations in the superfluid order parameter which are coupled to fluctuations in density. 
Therefore, we calculate the density response function as well as the pair response function. We use the GRPA which 
treats density and pairing fluctuations on an equal footing. In the context of ultracold Fermi superfluids, GRPA 
has successfully explained [33] Bragg scattering measurements [34] in the absence of an optical lattice. Various 
formulations of this method have been used to study the collective modes of Fermi superfluids in deep square/cubic 
optical lattices [24], [25|, [35|. Although GRPA is set up as a weak-coupling perturbation approach, it successfully 
captures the strong coupling limit as well. 

In this paper, we present GRPA results based on the Green's function approach developed by Cote and Griffin 
[36[. For the case of square/cubic lattice, this formulation [25| is in good quantitative agreement with a perturbation 
theory method used by two of us earlier [24]. For the honeycomb lattice case, we have explicitly checked that 
these two formulations are also in agreement. We briefly summarize the Cote- Griffin formalism in Appendix B. At 
strong coupling, we compare GRPA results with spin wave analysis of the appropriate pseudospin Hamiltonian. The 
formalism of the spin- wave analysis is summarized in Appendix C. 

A. Collective modes in the stationary ground state 

Figures EMU show the dynamic structure factor for density response and pairing response for different values of 
filling and flow. They are defined, respectively, as 

S d (q,w) = --\m X 7{q,m n ^Lo + i5), (15) 

7T 

S A (q,ou) = --Tmx%(q,iQ n ^LJ + i5). (16) 

7T 

Details are given in the Appendix B. The collective mode spectrum appears as a sharp peak in intensity of either 
dynamic structure factor, while the broad low intensity region is the two particle continuum. 

Since the GRPA equations of motion for the density and pair correlation functions have the same diagrammatic 
structure (see Eqs. (|B14|h (|B15|) , and Fig. [12]), they share the same poles and therefore the same collective mode 
spectrum. However, the spectral weight of the poles can be different for density and pairing correlation functions. 
The physical origin of collective modes can be understood from the relative weights for density and pairing response. 
To illustrate this, we discuss collective modes in the stationary ground state (Q = 0) in this section. 

The honeycomb lattice superfluid can have two undamped collective modes due to the presence of two sites in the 
unit cell, as shown in Fig. [6j 

(i) The lower mode incorporates the gapless Nambu-Goldstone mode corresponding to U(l) symmetry breaking. It is 
the usual 'Anderson-Bogoliubov' (AB) mode, predominantly composed of in-phase pairing fluctuations. This is seen 
from Fig. [6j where close to the T point, the AB mode shows a large intensity of the dynamic structure factor for 
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FIG. 6: (Color online) Calculated intensity of the (a) dynamic structure factor of density reponse function (S d (q,(jj)) and (b) 
that of pair response function (S A (q,uj)) in the ground state without flow (Q = 0) for U/t = 4, T = 0, and n = 0.9. The 
momentum q is along the contour displayed in the inset of Fig. The gapped Leggett mode (see text) has a large density 
component, while the gap less AB mode has a large pairing component. 



0.4 - 



0.2 




FIG. 7: (Color online) Comparison of the collective mode energy in strong coupling regime (U/t — 15) along the contour 
displayed in the inset. We set n = 0.5, T = 0, and Q = (0.2, 0). The GRPA result (solid line) is in good agreement with the 
spin- wave (S.W.) result (dashed line) for the strong-coupling pseudospin model. 



pairing response. 

(ii) The upper mode is a gapped 'optical' branch arising from density and/or phase fluctuations which are out of phase 
between two sublattices. This optical branch can be thought of as a 'Leggett mode' [37], in analogy with two-band 
superconductors such as MgE>2 [38j. When the interaction strength is not too large, as shown in Fig. [9](d), the Leggett 
mode may enter the two particle continuum and be strongly damped [l4j . 

Precisely at half-filling, the Hubbard model possesses a special SU(2) pseudospin symmetry [39|. Consequently, 
the superfluid ground state is degenerate with a CDW state which breaks sublattice symmetry. As a result, the 
Leggett mode which is composed of out-of-phase density fluctuations, becomes gapless. We are left with two Nambu- 
Goldstone modes, as the superfluid state breaks pseudospin SU(2) symmetry and not just U(l) gauge symmetry. 
Tuning away from half-filling, this degeneracy is weakly broken. The Leggett mode acquires a small gap as shown in 
Figs.[8](a)-(d)-(g). Close to half- filling, we deduce that the Leggett mode arises from low- lying CDW fluctuations - as 
shown in Fig. [6j this mode has much larger intensity for density response than for pairing. Far away from half-filling, 
CDW fluctuations are suppressed and both collective mode branches are predominantly composed of superfluid phase 
fluctuations. 

In addition to these collective modes, there exists a superfluid amplitude mode above the two particle continuum; it 
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FIG. 8: (Color online) Calculated intensity of the dynamic structure factor of density response function (S d (q,uj)) for different 
fillings and flow momenta, at U/t = 4 and T — 0. The momentum q is along the contour displayed in the inset of Fig. The 
Leggett mode goes down in energy with increasing flow and 'softens', signaling an instability. 



appears as a large peak at the edge of the continuum in Fig.[6](b). However, this amplitude mode is strongly damped 
by decay into pairs of single-particle excitations. Consequently, it may prove difficult to detect in experiments. We 
note that the amplitude mode very weakly couples with density fluctuations, if at all. As a consequence, the density 
correlation function does not show a peak at the edge of the continuum in Fig. [6] (a). 



B. Superflow instabilities 



As described in Sec. IIII A[ we impose superfluid flow on the system by forcing pairing at finite momentum. We 
use the GRPA method to extract the collective mode as a function of flow. As discussed earlier, GRPA also captures 
the strong coupling limit as shown in Fig. [7] by comparison with strong-coupling spin wave expansion. This gives us 
confidence that GRPA can be used to study superflow instabilities at all interaction strengths. Superfluidity can be 
sustained until some critical flow momentum which is set by one of three possible mechanisms [24j : 

(I) Depairing instability: the mean field superfluid order parameter vanishes. As discussed in Sec. IIII B) this instability 
occurs at the mean field level for weak to intermediate coupling strengths. 

(II) Landau instability: the collective mode energy hits zero and subsequently becomes negative [40|. The negative- 
energy modes are populated after the onset of the instability, leading to dissipation and the loss of superfluidity. As 
in the uniform superfluid Fermi gas [4lj |. this instability dominates in the strong-coupling BEC regime. 

(HI) Dynamical instability: the collective mode dispersion hits zero and subsequently becomes complex- valued. At 
the onset of this instability, one or more fluctuation modes start growing exponentially. Due to the rapid growth 
of fluctuations, this instability should be easy to observe in experiments. In fact, dynamical instabilities in bosonic 
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FIG. 9: (Color online) Calculated intensity of the dynamic structure factor of density response function (S d (q,uj)) for different 
fillings and flow momenta, at U/t = 4 and T = 0. The momentum q is along the contour displayed in the inset of Fig. [7] The 
Leggett mode undergoes a dynamical instability at the M point for n = 0.4. For n = 0.1, the usual Landau instability induced 
by long wavelength AB phonons occurs. 

condensates in a ID optical lattice have been successfully observed [42j . 

Figures [8] and [9] show the evolution of the collective mode spectrum as imposed flow is increased. As flow is increased, 
superfluidity is stable until one of the above instabilities occurs. Depending on filling and interaction strength, we 
can have different leading instabilities. Taking into account all the instabilities and their critical momenta, we map 
out a superflow 'stability phase diagram' in the intermediate coupling regime (U/t = 4) (Fig. [TO]) and strong coupling 
regime (U/t = 8) (Fig. [TT] left). We compare the strong coupling phase diagram obtained from GRPA with the one 
obtained from spin- wave analysis in the strong coupling regime (Fig. [TT] right) - see Appendix C for details. 

In the intermediate and strong coupling regimes, the collective mode spectrum always hits zero well before the two 
particle continuum does (see Figs. [8] and [9]). As a result, the critical flow momentum is always set by the collective 
mode as shown in the phase diagram in Fig. [10] The gapless superfluid phase and the depairing instability only exist 
at mean field level, and will play no role in experiments. They may become relevant to experiments at weak-coupling 
strengths. We discuss the instabilities seen in the stability phase diagram below: 

Dynamical instability at T: In the vicinity of half- filling (see Fig. |8]i), the Leggett mode has a small gap. This gap 
is a measure of the energy difference between the superfluid ground state and the low- lying sublattice-CDW phase. 
When a finite superflow is imposed, the kinetic energy of superflow overwhelms this excitation energy cost and the 
system switches to the CDW phase. This manifests itself as the dynamical instability of the Leggett mode at the T 
point. This is shown in Figs. [8] (d)~(f). As flow is increased, the Leggett mode goes down in energy and hits zero at 
the critical flow momentum. Subsequently, this mode energy becomes complex and sublattice-CDW fluctuations grow 
exponentially. We see this dynamical instability at the T point in the vicinity of and indeed precisely at half-filling 
(Fig.EI(a)~(c)). 

Landau instability: At low densities, since the system is similar to a continuum Fermi gas, the leading superflow 
instability is naturally a Landau instability of the AB mode. With imposed flow, the AB mode undergoes a 'Doppler 
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Filling 

FIG. 10: (Color online) Stability phase diagram for U/t = 4. The magnitude of the critical momentum is plotted as a function 
of flow imposed along the x-axis (T — K direction). In the light gray region, superfluidity is destroyed by collective mode 
instabilities despite the finite single-particle energy gap. The dark gray region indicates the gapless superfluid region (see text). 
Note that the depairing instability does not occur for 0.6 < n < 0.9 due to the robust superfluidity in the vicinity of the vH 
singularity at n = 3/4. 




Filling Filling 

FIG. 11: (Color online) Comparison of the stability phase diagram obtained by the GRPA (left) and spin- wave theory (right) 
in the strong-coupling regime U/t = 8. The magnitude of the critical momentum is plotted against flow imposed along the 
x-axis (r — K direction). 

shift' in the direction opposite to flow and becomes negative at long wavelengths (see panel (d)~(f) in Fig. [9j). This 
type of instability associated with negative energy collective modes has been observed in superfluid fermions in shallow 
optical lattices leading to the onset of dissipation (4lj . 

Dynamical instability at M: Surprisingly, at intermediate fillings (n c± 0.4), the leading superflow instability is 
a dynamical instability of the Leggett mode at the M point (see panels (a)~(c) in Fig. [9j). At low densities, this 
feature persists as a subleading instability occurring beyond the leading Landau instability (see Figs. [11]). At low 
densities, the unstable mode is dominated by phase fluctuations. At intermediate densities, it involves phase and 
density fluctuations in roughly equal measure. With our choice of superflow direction (Q along Y — K direction), this 
instability occurs at two M points which are symmetric with respect to the flow direction. This instability occurs even 
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when the flow is imposed along the Y — M direction; in which case, only one M point becomes unstable. The nature 
of the system beyond this instability is an exciting open question - we speculate that either local phase textures will 
form with local superfluid currents, or that the system will become chaotic. 

This M-point instability persists even in the strong coupling limit as shown in Fig. [TTJ This is in sharp contrast to 
the square lattice case where only Landau and sublattice-CDW instabilities occur in the strong coupling limit. In the 
strong coupling regime at low densities, we expect the physics to be dominated by phase fluctuations, captured in a 
Quantum Rotor (QR) model. In the square lattice case, the QR model is successful in predicting the first dynamical 
instability. This occurs at a flow momentum of tt/2 [43[ with the unstable wavevector at the BZ edge center. On 
the honeycomb lattice, the QR model predicts a dynamical instability at the flow momentum Q = (ir/y/3, 0) with 
the unstable wavevector at the M point. While it gets the wavevector of the instability right, it overestimates the 
critical flow momentum. Within GRPA and spin wave theory, we find that this instability occurs at Q ~ (27r/3\/3, 0). 
Since density fluctuations are not included in the QR model, they may suppress the critical flow momentum and give 
rise to the quantitative discrepancy. To clarify this issue, it may be useful to experimentally study the critical flow 
momentum for non-interaction bosons on a honeycomb lattice. For the square lattice case, such an experiment has 



Incommensurate dynamical instability: In Fig. EH the leading instability for intermediate fillings (0.5 < n < 0.7) at 
U/t = 4 is an incommensurate dynamical instability. This instability can be seen in Fig. [8] (g)~(i) with the unstable 
wavevector between the T and M points. Upon decreasing filling, the unstable wavevector moves towards the M 
point. Such an instability does not occur in the strong-coupling regime as shown in Fig. [TTJ This intermediate- 
coupling incommensurate instability has also been seen in square/cubic lattices 0, [25|. 

Figure [TTl right) shows the phase diagram in the strong-coupling limit from the spin wave theory. It exhibits 
qualitatively the same features as the strong-coupling phase diagram obtained by GRPA in Fig. Qjjleft): dynamical 
instability of the Leggett mode near the T point around half-filling, Landau instability of the AB mode at low density, 
and the dynamical instability of the Leggett mode at the M points for intermediate fillings. This demonstrates the 
consistency of the results obtained from GRPA and from spin wave analysis, serving as a check for the validity of 



In addition to the superflow instability discussed above, we point out some characteristic features of the spectra 
in Figs. [8]and [9j Close to half-filling, the Leggett mode is dominated by density fluctuations, while the AB mode is 
predominantly composed of superfluid phase fluctuations. This is clearly seen in Fig. [8] (e) and (h) from the small 
level repulsion between the two collective modes in the V — M region. At half-filling, as seen in Fig.[8](e) and (h), two 
branches of the collective mode overlap, indicating that their fluctuation components are orthogonal. On the other 
hand, at low fillings, since both collective modes are composed of superfluid phase fluctuations, they exhibit large 
level repulsion in Fig.[9](b) and (h) in the Y — M region. We also note that, at low densities, the Leggett mode enters 
the two particle continuum when a small flow is imposed as shown in Fig. [9j Consequently, the peaks become broad 
as they are damped by decay into single-particle excitations. However, the collective mode peaks still possess high 
intensity and appear as bright regions inside the continuum. 



We have studied the s- wave superfluid state on the honeycomb lattice inspired by recent cold atom experiments. The 
optical lattice setup of Ref. realizes a honeycomb geometry with tunable anisotropy. Starting from this configuration, 
we deduce that the symmetric lattice limit is optimal for the study of superfluidity. In addition, superfluidity is most 
robust at fillings close to n = 3/4, 5/4 fermions per site, due to a vH singularity in the non- interacting problem - this 
is the optimal density for observing superfluidity. Indeed, we find that the highest mean field occurs at the density 
corresponding to the vH singularity and in the symmetric lattice case. 

The critical velocity of superflow is a powerful tool to understand the low-lying excitations of a superfluid. To 
study this quantity, we first perform a mean field analysis of the superfluid state while allowing for imposed superflow. 
Experimentally, superflow can be imposed by using a 'running' optical lattice potential. On the honeycomb lattice, 
we find that superflow generically drives the system into a gapless superfluid state which is widest in the vicinity of 
the vH singularity at n = 3/4,5/4. Subsequently, the mean field order parameter vanishes in a first order depairing 
transition. 

We go beyond mean field theory by including density and pairing fluctuations using the GRPA formalism. At 
strong coupling, we compare our GRPA results with a spin-wave analysis of the appropriate pseudospin Hamiltonian. 
Although strictly justified only in the weak-coupling limit, GRPA works well even at strong coupling. We obtain 
sharp collective mode excitations with two branches - an Anderson-Bogoliubov mode and a Leggett mode. The 
former incorporates the Nambu-Goldstone mode from U(l) symmetry breaking. The latter is composed of out-of- 
phase fluctuations between the two sites of the honeycomb unit cell. At half-filling in the stationary superfluid, this 




GRPA. 
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Leggett mode becomes gapless at the T point, due to an extra pseudospin symmetry of the Hubbard model. This 
arises due to the sublattice-CDW state becoming degenerate with the superfluid. Away from half-filling, the Leggett 
mode acquires a gap reflecting the energy difference between the superfluid and the CDW states. 

With imposed superflow, the collective modes come down in energy due to a Doppler shift. An instability occurs 
when the collective mode softens - subsequently, the mode either becomes negative (Landau instability) or complex- 
valued (dynamical instability). Close to half- filling, the leading instability is a dynamical instability of the Leggett 
mode at the T point. This instability arises due to competition between the superfluid and sublattice-CDW phases. 
The presence of the low-lying CDW phase is a special feature of the one-band Hubbard model, arising from a delicate 
pseudospin rotational symmetry at half-filling. For future experiments, the presence of this instability close to half- 
filling can serve as a non-trivial test to verify if the one-band Hubbard model has indeed been realized. Naively, we 
may expect this instability to give rise to a stable 'supersolid' phase with coexisting superfluid and CDW orders. 
However, as in the square lattice case [23|, [24|, the honeycomb system is likely to go into a time-dependent non- 
equilibrium state. To probe this instability in experiments, it is best to focus on the short time window right after the 
critical momentum is reached. The exponentially growing CDW correlations can then be probed by rapidly ramping 
up the optical lattice to freeze the fermion density at each site. Subsequent imaging should give a snapshot clearly 
showing frozen CDW correlations. Even as we approach the critical flow momentum from below, we can use the same 
procedure to obtain snapshots which will show CDW fluctuations over some characteristic length scale. This length 
scale should diverge as we approach the critical flow momentum. 

At low densities, we find a Landau instability of the AB mode. We also find an unexpected dynamical instability 
at the M point for intermediate densities. Surprisingly, this instability persists even at strong coupling. The nature of 
the system beyond this instability is an interesting open question. We present stability phase diagrams which give the 
critical flow momentum and the associated instability mechanism as a function of density and interaction strength. 

Our calculations are also relevant to recent studies of the repulsive Hubbard model on the honeycomb lattice. At 
half- filling, quantum Monte Carlo simulations (45| have suggested a gapped spin liquid state at intermediate U/t 
between semi- metal and the Neel phases. There have been many attempts to understand this intermediate phase 
and the nature of phase transitions into and out of it. Precisely at half-filling, the repulsive and attractive Hubbard 
models are related by a sublattice-dependent particle hole transformation. In fact, the attractive model at any filling 
is equivalent to the repulsive problem at half- filling, but with a magnetic field [46[. Our GRPA approach does 
not capture a strongly correlated spin-liquid phase as suggested by Quantum Monte Carlo simulations. However, 
approaching the spin liquid from the Neel side, the Neel-spin liquid quantum phase transition may arise from melting 
of Neel order due to collective mode excitations. Such a mechanism has been discussed in Refl471 in the context of 
the Bi3Mn40i2NO(3), which is a spin-3/2 honeycomb lattice antiferromagnet. The question can be phrased within 
the attractive Hubbard model - can collective modes obtained from GRPA melt superfluid order? If so, what is the 
nature of the resulting state? An exciting possibility is to obtain a 'pairing liquid' state which is the analog of the spin 
liquid proposed for the repulsive model [49j. Such a pairing liquid can be thought of as a quantum pseudogap phase 
with preformed Cooper pairs. Recently, a finite-temperature pseudogap phase has been observed in a trapped Fermi 
gas near the unitary limit [5(J However, a zero-temperature analog has not been explored. This is an interesting 
direction for future research. 
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Appendix A: Mean field theory equations 

The Hamiltonian of Eq. ([7]) can be solved by a Bogoliubov transformation as given in Eq. (j9j). Denoting the elements 
of the transformation matrix T by (T%j = T^, the gap equation (j6j) gives 

U 



A Q ,a = ^ E [T&TnfiE+ik)) + T* 2 T 12 (l - /(£+(-*)) 

k 

+T 2 * 3 T 13 /(^(fe)) + T 2 * 4 T 14 (l - f(EQ(-k)) 



A Q,b = 



U_ 

M 



[T^fiE+ik)) + Ti 2 T 32 (l - f(E+(-k)) 

k 

+T: 3 T 33 f(EQ(k))+T; 4 T 34 (l - /(££(-fe)))' , 

where (a£ r(7 ttfe )T)(T ) = f(Eq(k)) (f(E) = 1/je^ + 1} is the Fermi distribution function). 
The average filling per site n is 

n = ^EE< nfc ^> 

k 

= jf E [(l T ll| 2 - l^ll 2 + l T 3l! 2 - l T 4l| 2 ) f(E+(k)) 

k 

+ (~\T 12 \ 2 + |T 22 | 2 - |T 32 | 2 + |T 42 | 2 ) f(E+(-k)) 
+ (|T 13 | 2 - |T 23 | 2 + |T 33 | 2 - |T 43 | 2 ) /(^(fc)) 



(Al) 



(A2) 



(-lH- 



|T 34 | 2 + |T 44 | 2 )/(i^(-fc)) 



|Ti 2 | 2 + |T 14 | 2 + |T 21 | 2 + |T 23 | 2 + |T 32 | 2 + |T 34 | 2 + |T 41 | 2 + |T 43 | 2 1 



(A3) 

With imposed superflow (Q ^ 0), these equations have to be solved numerically For the stationary superfluid, 
however, setting Aq iCL = Aq ^ = Ao, the T matrix can be obtained analytically 

\ 



Tq=0 = 71 



+ 



V 



,+ 



-i<Pk 



k 



(A4) 



) 



where e^ fe = Jk/\jk\- The matrix elements u k and v k in Eq. (|A4|) have the standard form as 



E±{k) 



(A5) 



The gap equation ()A1|) and (|A2|I reduce to 



?7 



A ° = ]vEE«[l-2/(£ T (fe))] 



ua ^2 ^2 tan ^ 



TV 



k T=± 



2E T (k) 



Thus, we obtain 



The number equation (|A3|) becomes 



/7 N ^ ^ 



tanh 



/3£T(fc) 



2E T (k) 



(A6) 



(A7) 



= ^ E E KK) 2 - Mt) 2 }/(^(*0) + K) 2 ] 



r=± 



i-IyyJ 



tanh 



/3E T (k) 



(A8) 
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The above equations (|A7|) and (|A8|) have been derived in Ref. [l4j . 

Appendix B: GRPA formalism of Cote and Griffin 

In this Appendix, we give a brief summary of the GRPA formalism developed by Cote and Griffin [36]. In this 
formalism, we introduce fictitious external fields to calculate various response functions by taking functional derivatives 
of the single-particle Green's function with respect to the external fields. The collective mode spectrum of the system 
is given by the poles of the response functions. As emphasized in Ref. [36], this technique has the advantage that 
the response functions are uniquely determined by the functional differentiation technique to be consistent with the 
self-energy approximation, and therefore gives a gapless Nambu-Goldstone mode in the superfluid phase. 

The Hamiltonian with time-dependent external fields coupled to density and pair operators is 

K{t) =H + Y j Pi(r)ni + (QiWciW + Qi(r)cl4 t ) , (Bl) 

i i 

where r is imaginary time. Hereafter, we adopt the notation 1 = (ji, n) = (Zi, ti), where j\ denotes the lattice site, 
while l\ and v\ denote the unit cell and sublattice, respectively. It is convenient to introduce the Nambu representation 
[48| , in which a single-particle Green's function is given in a 2 x 2 matrix form, 



(B2) 



G(l,2) = -(T*(l)* t (2)) 

/ -e- i ^-( ,, 'i- ri 2)(Tc t (l)c|(2)) -e-^-( r! i +,,i 2)(Tc t (l)c4.(2)) \ 
^ -e i ?-(^i+n 2 )(Tc|(l)c|(2)) e i T'(n 1 -n 2 )^ Tc; ( 2 ) c t(i)) y 

The field operator is a spinor consisting of the creation and annihilation operators 

= ^ ( 4(1) ) ' * t(1) = O^ 1 )' C + W) "ft" ( B3 ) 

The unitary matrix 7^ is given by 



We note that the Green's function in Eq. (|B2[) is defined in the co-moving frame of the superfluid. The Green's 
function in the laboratory frame is obtained by the unitary transformation 7^(3(1, 2)7^. 
The Green's function G(l,2) satisfies the Dyson equation 

\G{\M~ X = [G^M)]- 1 - £(1,2) - W(l,2), (B5) 
where G°(l,2) is the non-interacting Green's function. The external field matrix W(l,2) is defined as 

W(l, 2) - 6(1 - 2) ( _ e Z%, {1) - e_ : Q ; ( i ; ) g(1) ) ■ (B6) 
The self-energy within the Hartree-Fock-Gor'kov (HFG) approximation is given by 

E(l,2) = *(l-2)(:^ (B7) 
We introduce the three-point density correlation function as 

i(l,2,3 )S §^>, ,B8, 

where G(l,2) = 73(5(1,2) with T3 being the usual Pauli matrix. In the equal space-time limit, it reduces to 

L(1,3) = L(1,1 .-J)-^ e «}.r ll(Tftn t(i) jn (3)) -(T5 ni (l)6n(3)) J' 



17 








FIG. 12: The diagrammatic representation of the GRPA equations (|B14|) and (|B15[) . The three point correlation function A 
is either the density (X) or pair (M) correlation function. The wavy line represents the attractive interaction —U. A involves 
the bubble diagrams, while the irreducible correlation function A involves the ladder diagrams. 



The density response function is thus obtained from the diagonal elements of Eq. (|B9|) as 

Xd(1 ' 3) = jp^ = -(T5n(l)Sn(3)) 
= L 11 (l,3) + L 22 (l,3). 
We also introduce the three-point pair correlation function as 



M(l,2,3) 



SQ(3) ■ 



In the equal space-time limit, it reduces to 



3(T<5n t (l)(5mt(3)) - e -"'-( r 'i- p 's) (TSm(l)Sm^(3)) 



e«J-( r 'i-' , «s)(T<5mt(l)(ymt(3)> -e^'s (TSn^Sm* (3)) 



M(l,3) = Af(l,l+3) = 
The pair response function is obtained from the off-diagonal element of Eq. (jB12[) as 



(BIO) 



(Bll) 



(B12) 



Xa(1,3) 



5Q(3) 

iQ < ri i- ri ^{T5m(l)5m\3)) = M 12 (l,3). 



(B13) 



Differentiating (5(1, 2) with respect to the external fields and using the HFG approximation, one obtains the GRPA 
equation of motion for the three-point response functions as 



4(1,2,5) = 1(1, 2, 5) -U J d3L(l,2,3)[An(3,5) + A 22 (3,5)], 
1(1,2,5) = i°(l,2,5) + Z7 / d3 G(l, 3)1(3, 5)G(3, 2), 



where A is either L or M. The lowest-order correlation functions are given by 



L°(l,2,5) 
M°(l,2,5) 



G(1,5)G(5,2), 
G(l,5)(|] I ) 



(B14) 
(B15) 

(B16) 
(B17) 



A(l,2,3) is called the irreducible correlation function. The form of Eqs. (|B14|) and (|B15|) exposes the diagrammatic 
structure of the GRPA, as shown in Fig. [I2j where Eqs. (|B14[) and (|B15[) represent the summation of the ladder 
diagrams and bubble diagrams, respectively. 
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The GRPA equations in momentum space are 

A Ul ^(q) = A u ^(q)-U^2L u ^(q)A u ^(q), (B18) 

A™(q) = A^W + ^J^Yt^iP^^MG^ip-q), (B19) 

where q = (q,iQ n ) and p = (p,iuj n ) (uo n and ^ n are Fermi and Bose Matsubara frequencies, respectively). Here, we 
use the Fourier transform 

G VlV2 (l,2) = -j- J2 G UlU2 (p)e^ p < ri i- r ^-^^- T ^\ (B20) 

A" lV2 (l,2) = -L ^ i^^^)^(q-(rz 1 -n 2 )-^(r 1 -r 2 )) > ( B21 ) 

Using Eqs. (|B1Q|) and (jB13|h the density and pair response functions can be evaluated from the correlation functions 
L and M that are obtained by solving Eqs. (|B18[) and (|B19j) . 

Appendix C: Spin-wave analysis in the strong coupling limit 

In the limit of large U/t, states with singly occupied sites can be ignored leading to an effective spin- 1/2 Heisenberg 
model [24l ]. The effective Hamiltonian can be written as 

H=J S ha -S t+Stb -(BS) S tv ( C1 ) 

i, 5=0,a,2 ,ai+a2 i,v=a,b 

The exchange interaction J = At 2 /U arises from superexchange. In this spin mapping, the X and Y directions of 
spin correspond to the local superfluid order parameter. The density of fermions maps onto the z-polarization, with 
half-filling mapping to zero z-polarization. The chemical potential in the fermion problem gives rise to a magnetic 
field tuning polarization. We take the magnetic field to be proportional to S, so that all terms in the Hamiltonian 
scale as S 2 . 

The uniform superfluid state at half-filling maps onto a Neel antiferromagnetic order with spins lying in the XY 
plane. This state is symmetric under rotation about the Z axis, which corresponds to the usual U(l) gauge symmetry. 
Superflow is imposed by a phase gradient, which corresponds to twisting the local Neel vector at by an angle Q • 
leading to a spiral state with spins lying in the XY plane. Away from half-filling, a non-zero magnetic field is 
introduced. This forces the spins to cant uniformly towards the Z axis to gain Zeeman energy. We characterize the 
flowing superfluid as 

S;, a = S'(cosxcos(Q.r i ),cosxsin(Q.r i ),sinx) 
Si, 6 = S(- cosxcos(Q.r i + 0), - cosxsin(Q.r i + 0),sinx) (C2) 

We allow the A and B sublattices to have a relative phase difference (j). The Neel state (stationary superfluid) is 
recovered by setting Q = <fi = 0. For simplicity, we denote the phase angles by oti = Q • and ft = Q • + <p. We 
now perform a local spin rotation and define new spin operators which we denote as T iyl/ = Wi jU Si jU , by 

Tia\ ( sm X — cosx\ / cos(c^) sin(c^) 
Tf a = 1 -sin(a<) cos(a,) 

TtJ \ cos X sinx / V 01 




1 i,b 



Lb 



— sin x — cos x 

1 
cos x — sin x 



cos(ft) sin(ft) 0\ /Sf ih 
-sin(A) cos(ft) ) ( Sl h 











(C3) 



J Lb 



With this local spin rotation, the flowing superfluid state transforms into a uniform antiferromagnet with all spins 
pointing towards ±Z. We absorb the spin rotation operators into the Hamiltonian. The exchange coupling then takes 
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the form 



Si a ' S 



3,b 




S Y S ( 



-s 2x c 2 ai 



-Pj)/2 C X S Q 



ati—/3j 



S ^X C la x -P 3 )/2 
-C x S ai -f3 3 







II 









-Pj 4 u ^-ft- -x 

where j = i + 5. We denote sin (9 = Se and cos (9 = Cq. The Zeeman term in the Hamiltonian becomes 

- BS (St a + S? 6 ) = —BS {-C x {T* a + T? b } + S x {T? a - T? b }) . 
We now introduce Holstein Primakov bosons as follows: 



(C4) 



i,a 
i,a 



(a.i + 4), 

f (*)(«* -<4). 



1 i,b 



'§ (t) Qh - b\), 

b]h - S. 



The bosonic operators a, and 6, represent fluctuations about the spiral state. S is the spin length, which we will 
ultimately set to be S = 1/2. The Hamiltonian can now be expanded in powers of S. The leading terms, of order S 2 , 
give the energy in the classical limit. 



Eci = 2NS 2 



-JC 



2lQ e 



7-Qe" 



JS^€q — 2BS X 



(C5) 



where TV is the total number of sites. We have denoted 7k = 1 + e zk a2 + e «k-(ai+a 2 )^ Terms proportional to S 3 ^ 2 are 
linear in boson operators. The coefficient of the single boson operators is complex- valued - we set this quantity to be 
zero. Setting the real part to be zero, we get 



smx 



B 



J(3 



2 



) 



(C6) 



Alternatively, this equation be obtained by treating x as a variational parameter and demanding d x Eci = 0. Setting 
the imaginary part of the coefficient to zero gives 

sin(Q 2 ) +sin(Qi + Q 2 ) 

tan (7) = — — — — — , (C7) 

Y l + cos(Q 2 ) + cos(Qi + Q 2 ) V ; 

where Qi and Q 2 are the components of Q along the primitive vectors a\ and a 2 shown in Fig. [U This fixes the 
value of <j). Again, this equation can be also obtained variationally by demanding d^Eci = 0. 

The next terms in the Holstein-Primakov expansion are of order S. These capture the quantum spin wave fluctua- 
tions. These terms are given by 



H 



2NS 



3JS 2 X -JC 2 X ^- 



-7-qe 



BS Y 



where 



b k 


( 




V bl k J 


V 



X k + F k Z k 
Uk Z k X k — Y k 

u k 



(C8) 



As the Hamiltonian matrix is Hermitian, we have only filled in the upper triangle. The entries are given by 



X k 



J(-3S* + C 



X 



itfr 



7-Qf 



-)+BS x , 



C x 7k - (1 + S x ) 



-JS 7 
2 * 



Tk+qe 6 



7k-Qe 



-i(f) 



7k' 



2 

7k+qe i 



7k-Qe" 



(C9) 
(CIO) 
(Cll) 
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We obtain the spin wave spectrum by performing a bosonic Bogoliubov transformation on this Hamiltonian matrix. 
For a given filling and flow momentum, we can thus obtain the collective mode excitations in the strong coupling 
limit. 
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